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Q ! ABSTRACT 

ryj t In the present paper, we improve the "Extended Secondary Infall Model" (ESIM) of Williams 
C3 . et al. (2004) to obtain further insights on the cusp/core problem. The model takes into account 
the effect of ordered and random angular momentum, dynamical friction and baryon adiabatic 
, contraction in order to obtain a secondary infall model more close to the collapse reality. The 
• model is applied to structures on galactic scales (normal and dwarf spiral galaxies) and on 
' cluster of galaxies scales. The results obtained suggest that angular momentum and dynamical 
, friction are able, on galactic scales, to overcome the competing effect of adiabatic contraction 
t-H ■ eliminating the cusp. The NFW profile can be reobtained, in our model only if the system is 
' constituted just by dark matter and the magnitude of angular momentum and dynamical friction 
i are reduced with respect to the values predicted by the model itself. The rotation curves of four 
f"^ 1 LSB galaxies from de Blok & Bosma (2002) are compared to the rotation curves obtained by 
^\ t the model in the present paper obtaining a good fit to the observational data. On scales smaller 
. than ~ 10 11 /i _1 Mq the slope a ~ and on cluster scales we observe a similar evolution of the 
^ ' dark matter density profile but in this case the density profile slope flattens to a ~ 0.6 for a 
• I \ cluster of ~ lO u ft _1 M . The total mass profile, differently from that of dark matter, shows a 
■ central cusp well fitted by a NFW model. 

b ! 

. . i Key words, cosmology-theory-large scale structure of Universe-galaxies-formation 



1. Introduction 

The structure of dark matter haloes is of fundamental importance in the study of the 
formation and evolution of galaxies and clusters of galaxies. At the simplest level, dark 
matter halos form when, in the early Universe, the matter within (and surrounding) 
an overdense region suffers gravitational retardation, decouples from the Hubble flow, 
collapses, and in due course, virializes. From the theoretical point of view, the structure 
of dark matter haloes can be studied both analytically and numerically. A great part of 
the analytical work done so far is based on the secondary infall model (SIM) introduced 
by Gunn & Gott (1972), Gott (1975) and Gunn (1977). Calculations based on this model 
predict that the density profile of the virialized halo should scale as p oc r~ 9 / 4 . Self-similar 
solutions were found by Fillmore & Goldreich (1984) and Bertschinger (1985), while 
Hoffman & Shaham (1985) (hereafter HS) studied density profiles around density peaks. 
More recently modifications of the self-similar collapse model to include more realistic 
dynamics of the growth process have been proposed (e.g. Avila-Reese et al. 1998; Nusser 
& Sheth 1999; Henriksen & Widrow 1999; Subramanian et al. 2000; Del Popolo et al. 2000 
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(hereafter DP2000)). Rydcn & Gunn (1987) (hereafter RG87) were the first to relax the 
assumption of purely radial self-similar collapse by including non-radial motions arising 
from secondary perturbations in the halo. Numerous authors have emphasized the effect 
of an isotropic velocity dispersion (thus of non-radial motion) in the core of collisionless 
haloes. One common result of the previous studies (RG87; White & Zaritski 1992; Avila- 
Reese et al. 1998; Hiotelis 2002; Nusser 2001; Le Delliou & Enriksen 2003; Ascasibar 
et al. 2004; Williams et al. 2004) is that larger amount of angular momentum leads to 
shallower final density profiles in the inner region of the halo. Moreover, baryons have 
been invoked both to shallow (El-Zant et al. 2001 (hereafter EZ01), 2003; Romano-Diaz 
et al. 2008) and to steepen (Blumenthal et al. 1986) the dark matter profile. 

As previously reported, the structure of dark matter haloes can also be studied 
through numerical simulations. Quinn, Salmon & Zurek (1986) pioneered the use of N- 
body simulations to study halo formation and confirmed HS results. More recent studies, 
(Dubinski & Carlberg 1991, Lemson 1995, Cole & Lacey 1996, Navarro et al. 1996, 1997 
(NFW), Moore et al. 1998, Jing & Suto 2000, Klypin et al. 2001, Bullock et al. 2001, 
Power et al. 2003 and Navarro et al. 2004) found that although the spherically-averaged 
density profiles of the N-body dark matter halos are similar, regardless of the mass of the 
halo or the cosmological model, their profiles are significantly different from the single 
power laws predicted by the theoretical studies. The N-body profiles are characterized by 
an r -3 decline at large radii and a cuspy profile of the form p(r) oc r~ a , where a < 2 near 
the center. The actual value of the inner density slope a is a matter of some controversy, 
with NFW suggesting a = 1, but with Moore et al. (1998), Ghigna et al. (2000) and 
Fukushige & Makino (2001) arguing for a = 1.5, while Jing & Suto (2000) and Klypin et 
al. (2001) claimed that the actual value of a may depend on halo mass, merger history, 
and substructure. Power et al. (2003) pointed out that the logarithmic slope becomes 
increasingly shallow inwards, with little sign of approaching an asymptotic value at the 
resolved radii. In that case, the precise value of a, at a given cut-off scale, would not be 
particularly meaningful. This result has been later confirmed by Hayashi et al. (2003) 
and Fukushige et al. (2004), and it is predicted by several analytical models (e.g., Taylor 
& Navarro (2001), Hocft et al. 2003). Finally, Navarro et al. (2004) proposed a new fitting 
formula having a logarithmic slope that decreases inward more gradually than the NFW 
profile. 

While numerical simulations universally produce a cuspy density profile, observed 
rotation curves of dwarf spiral and LSB galaxies seem to indicate that the shape of the 
density profile at small scales is significantly shallower than what is found in numerical 
simulations (Flores & Primak 1994; Moore 1994; Burkert 1995; Kravtsov et al. 1998; 
Salucci & Burkert 2000; Borricllo & Salucci 2001; de Blok et al. 2001; de Blok & Bosma 
2002; Marchcsini et al. 2002; de Blok 2003; de Blok, Bosma & McGaugh 2003). It seems 
that the data generally favor logarithmic density slopes close to 0.2 (de Blok 2003; de 
Blok, Bosma & McGaugh 2003; Spekkens et al. 2005). 

On cluster scales, X-ray analyses have led to wide ranging results, from a = 0.6 
(Ettori et al. 2002) to a = 1.2 (Lewis et al. 2003) or even a = 1.9 (Arabadjis et al. 
2002). Ricotti's (2003) N-body simulations suggest that density profile of DM haloes is 
not universal (in agreement with Jing & Suto 2000; Subramanian et al. 2000; Simon et 
al. 2003b; Cen et al. 2004; Ricotti & Wilkinson 2004; Ricotti et al. 2007), presenting 
shallower cores in dwarf galaxies and steeper cores in clusters. 

The discrepancy between simulations and observations has been often signaled as a 
genuine crisis of the CDM scenario and has become known as the "cusp/core" problem. 
Since LSB galaxies are thought to be ideal for the comparison with theory, as their 
dynamics are dominated by dark matter with little contribution from baryons (Bothun et 
al. 1997), the discrepancy with simulations is particularly troublesome. The significance 
of this disagreement, though, remains controversial and different solutions have been 
proposed. A number of authors attribute the problem to a real failure of the CDM 
model, or to that of simulations (de Blok et al. 2001a; de Blok, McGaugh, & Rubin 
2001b; Borricllo & Salucci 2001; de Blok, Bosma, & McGaugh 2003). This has led to 
suggestions that dark matter properties may deviate from standard CDM and several 
alternatives have been suggested, such as warm (Colin et al. 2000; Sommer-Larsen & 
Dolgov 2001), repulsive (Goodman 2000) fluid (Peebles 2000), fuzzy (Hu et al. 2000), 
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decaying (Cen 2001), annihilating (Kaplinghat et al. 2000), or self-interacting (Spergel 
& Steinhardt 2000; Yoshida et al. 2000; Dave et al. 2001) dark matter. Others argue 
that the inconsistency may reflect the finite resolution of the observations that has not 
been properly accounted for in the analysis of the HI rotation curves (van den Bosch et 
al. 2000; van den Bosch & Swaters 2001; Rhee et al. 2004). Alternatively, it has been 
suggested that stellar feedback from the first generation of stars formed in galaxies was 
so efficient that the remaining gas was expelled on a timescale comparable to, or less 
than, the local dynamical timescale. The dark matter subsequently adjusted to form 
an approximately constant density core (e.g., Gelato & Sommer-Larsen 1999). This is 
however unlikely to affect cluster cusps. 

Simulations, observations and semi-analytic models agree on outer parts of the haloes' 
structure the disagreement in the inner regions could be connected to limits in numerical 
simulations (de Blok 2003; Taylor et al. 2004) or to the fact that dissipationless N-body 
simulations do not take into account the effects of baryons on dark matter evolution. 
Interestingly, the amount of central substructure seen in the semi-analytic haloes is con- 
sistent with the amount of substructure inferred from strong lensing experiments (Taylor 
et al. 2004). Thus the semi-analytic haloes might provide a more accurate picture of 
the spatial distribution of substructure around galaxies even if an analytical method, no 
matter how sophisticated, will never be able to capture the full extent of complexity of 
a non-linear process. 

Nevertheless the very large amount of work carried out by many researchers to date 
using N-body simulations has met with limited success in elucidating the physics of halo 
formation. The reason is due to the fact that the point of force of numerical simulations 
(namely to capture the full extent of complexity of a non-linear process) is also their 
weakness: numerical simulations yield little physical insight beyond empirical findings 
precisely because they are so rich in dynamical processes, which are hard to disentan- 
gle and interpret in terms of underlying physics. Analytical and semi-analytical models 
are much more flexible than N-body simulations (see Williams et al. 2004). So even 
if analytical models like SIM treat collapse and virialization of halos that are spheri- 
cally symmetric, that have suffered no major mergers, and that have suffered quiescent 
accretion, they are worth investigating. One of the most used semi-analytical models 
is the SIM whose most often questioned assumptions are the spherical symmetry and 
the absence of peculiar velocities (non-radial motions): in the "real" collapse, accretion 
does not happen in spherical shells but by aggregation of subclumps of matter which 
have already collapsed; a large fraction of observed clusters of galaxies exhibit signifi- 
cant substructure (Kriessler et al. 1995). Motions are not purely radial, especially when 
the perturbation detaches from the general expansion. Nevertheless the SIM gives good 
results in describing the formation of dark matter haloes, because in energy space the 
collapse is ordered and gentle, differently from the chaotic collapse that is seen in N- 
body simulations (Zaroubi, Nairn & Hoffman 1996). This is confirmed in other studies 
(Toth & Ostriker 1992; Huss, Jain & Stcinmctz 1999a,b; Moore et al. 1999). We should 
also add that analytical and semi-analytical models have some advantages over N-body 
simulations: a) they are flexible (one can study the effects of physical processes one at a 
time); b) one can incorporate many physical effects at least in a schematic manner; c) 
they are computationally efficient (it takes about 10 s to compute the density profile of 
a given object at a given epoch on a desktop PC, Ascasibar et al. 2007). 

In this paper, we shall present an analytical model for haloes formation based on SIM. 
The model is an extension of the Williams et al. (2004) model to take account of ordered 
angular momentum, dynamical friction and adiabatic contraction. The model derives the 
initial shape of proto-halos from the fluctuation spectrum at high redshifts, and halos are 
endowed with secondary perturbations, which impart random motions to halo particles, 
as in Williams et al. (2004), and ordered angular momentum. The statistical properties 
of the secondary perturbations are derived from the same fluctuation power spectrum, 
and therefore their effects on particles random velocities are treated self-consistently. 
The ordered angular momentum acquired by the proto-structures is obtained from the 
tidal interaction of the proto-structure with the neighboring ones and using the theory of 
random fields. Dynamical friction is calculated according to Kandrup's (1980) approach. 
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Adiabatic contraction is taken into account using the Blumenthal et al. (1986) model as 
modified by Gnedin et al. (2004). 

The plan of the paper is the following: in Section 2, we introduce the model that shall 
be used to calculate the density profiles. In Section 3, we show the results of the model, 
while Section 4 is devoted to conclusions. 



2. Model 

As shown by the spherical collapse modeQ of Gunn & Gott (1972), a bound mass shell 
having initial comoving radius Xi will expand to a maximum radius x rn (named apap- 
sis or turn-around radius Xt a ). As successive shells expand to their maximum radius 
they acquire angular momentum and then fall in on orbits determined by the angular 
momentum. Dissipative physics and the process of violent relaxation will eventually in- 
tervene and convert the kinetic energy of collapse into random motions (virialization) . 
Following RG87 and R88a,b aproach, we restrict ourselves to halos that form primarily 
via nearly-smooth accretion of matter, treating collapse and virialization of halos that are 
spherically symmetric, that have suffered no major mergers, and that have experienced 
only quiescent accretion of somewhat lumpy material. The most important feature of the 
original RG87 method is that the dynamical evolution is carried out while conserving 
angular and radial momenta of individual halo shells. The initial shape of the proto- 
halos is derived from the fluctuation spectrum at high redshifts, and halos are endowed 
with secondary perturbations, which impart random motions to halo particles. The sta- 
tistical properties of the secondary perturbations are derived from the same fluctuation 
power spectrum, and therefore their effects on particles random velocities are treated 
self-consistently. Ordered angular momentum is also treated self-consistently as shown 
in R88a. 

It is possible to divide the halo into four regions starting from the outside: in 
Region (1), furthest from the center of the initial density peak, dark matter particles 
are beginning to feel the gravitational tug of the central peak and are just starting 
to fall behind the Hubble flow. Further, in Region (2), the particles are starting to 
decouple from the Hubble flow and are about to begin collapse. In Region (3) the central 
density peak dominates the motion of particles; this is the region of infall and shell 
crossing. Finally, in Region (4), the central part of the density peak, virialization is 
taking place, or has already been reached. In Region (4), virialization is well underway 
or is complete. Many analytical solutions in regions (l)-(3) ignore non-radial motions of 
particles; all orbits are assumed to be purely radial. The calculations carried out in RG87 
demonstrate that non-radial motions are very important in determining the outcome in 
this region. Nonradial motions give particles angular momentum, which prevents them 
from penetrating to the very center of the halo (confirmed by e.g., Avila- Reese, Firmani 
& Hernandez (1998), Subramanian et al. (2000), and Hozumi, Burkert & Fujiwara (2000). 

In most popular cosmological scenarios the density field soon after recombination 
can be represented by a Gaussian random field. High density contrast peaks in the field 
will eventually achieve overdensities of order 1 and enter a non-linear stage of evolution. 
These peaks will then collapse to form bound structures. We start with one of these 
peaks, and, for simplicity take it to be spherically symmetric. The peak is divided into 
a very small central core and many spherically symmetric concentric mass shells, each 
labeled by its initial comoving distance from the center, x. 

The first step in obtaining the density profile, is to calculate the initial density profile 
arising from a primordial fluctuation, 8i(xi). A well known result is the expression for the 



1 A slightly overdense sphere, embedded in the Universe, is a useful non-linear model, as it 
behaves exactly as a closed sub-universe because of Birkhoff 's theorem. The sphere is divided 
into spherical "shells" . A spherical "shell" may be defined as the set of particles at a given radius 
that are all at the same phase in their orbits (see Le Delliou & Henriksen 2003). 
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radial density profile of a fluctuation centered on a primordial peak of arbitrary height 



(5(x)} = S (x) 



v£(x) J?(^7>7) 



m i/2 7 (1_ 7 2) 



2 



m~ 1/2 (i) 



(BBKS; RG87), where x is the comoving separation, v = 5(0) /(^ is the height of a 
density peak, £(r) is the two-point correlation function: 

£(x) = -^-[ P(k)ksm(kx)dk (2) 
2tt z x J 

7 and i?* are two spectral parameters (see BBKS) and $(^7,7) is given in BBKS. 

The CDM spectrum used in this paper is that of BBKS (equation (G3)), with transfer 
function: 

T(k) = [ln(1 2 + 34 2 g 34g)] • [1 + 3.89q + (W.lqf + (5A6qf + (6.71) 4 ]- 1/4 (3) 

where q = iix ^f M / pc -i ■ Here 9 = p cl -/ (1.68p 7 ) represents the ratio of the energy density in 
relativistic particles to that in photons (8 = 1 corresponds to photons and three flavors 
of relativistic neutrinos). The spectrum is connected to the transfer function through the 
equation: 

P(k) = PcDMe-^ k2R2 f (4) 

where Rf is the smoothing (filtering) scale and Pcdm is given by: 

P C DM=AkT 2 (k) (5) 

where A is the normalization constant. We normalized the spectrum by imposing that 
the mass variance of the density field 

1 f°° 

a 2 (M) = — / dkk 2 P(k)W 2 (kR) (6) 
convolved with the top hat window 



3 

W(kR) = — — (sin kR - kR cos kR) (7) 
(kR) 6 

of radius 8 hr 1 Mpc -1 is as = 0.76 (Romano-Diaz et al. 2008). Throughout the paper we 
adopt a ACDM cosmology with WMAP3 parameters, Q m = 1 - Oa = 0.24, Q\ = 0.76, 
fib = 0.043 and h = 0.73, where h is the Hubble constant in units of 100 km s _1 Afpc -1 . 
The density profile of an initial, pre-collapse halo, at early times, Si(x), is related to 5q(x) 
by the linear growth factor D(z) (Peebles 1980) by: 

Si(x) = 8 (x)/D( Zi ) (8) 

Let us divide the main smooth spherically symmetric halo, 5(x), into many concentric 
mass shells. Each shell is uniquely labeled by x, its initial comoving radius. The halo 
represents an upward departure from average background density, i.e. density interior 
to any shell is greater than critical at all times. Each shells evolution is divided into 
two stages. Initially a shell expands with the Hubble flow, but with a slight deceleration 
arising from the central mass concentration. Eventually the shells outward radial velocity 
decreases to zero, after which the shell collapses, by some distance, back towards the 
center of the halo. The dividing moment is called the turn-around, and at any given time 
corresponds to the line dividing Regions (1) and (2). In a halo with a declining density 
profile the turn-around happens at progressively later cosmic times for shells at greater 
distances from the center. 



2 a is the mass variance filtered on a scale Rf 
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In such cases the time evolution of a shell until it reaches turn-around is given by a 
set of parametric equations (Gunn & Gott 1972, Peebles 1980), 

r(0) = lx5o 1 (l-a»0), (9) 

t{0) = ^t 6o 3/2 (9-sm6), (10) 

where r is particle's proper radius, t is cosmic time, and So(x) is the initial average 
fractional density excess inside the shell (Eq. [13] of RG87), 

= \ I S Q (y)y 2 dy, (11) 
x Jo 

and to is the present time. The mass of a shell, and the mass within a shell are constant; 
these relations, together with Eqs. [9] and [10] can be used to compute the fractional 
overdensity for any shell, parameterized by its So /So, at a cosmic time corresponding 
to 9, (Eq. [23] of RG87): 



9(6»- sin 0) 2 



So 



xta\ j.i " ' i 1 i p 1 ■» i 3 sin 0(0 - sin 0) i \ 

<5W + 1= 2(l-cos0)3l v 1 + 3 [ 1 -^J[ 1 - 2(l-cos0)2 \) ' (12) 

This can be used to construct the density run of a halo at a fixed cosmic time t/to- 

The mean density distribution about a peak is spherically symmetric. However, the 
initial density peak will in general have a triaxial shape, leading to non-spherical collapse. 
As shown by R88a, since the quadrupole moment of the protogalaxy is largely due to the 
outermost shells, where triaxiality is insignificant, it is justifiable to ignore the triaxility 
in computations. In other terms, the dynamics of the halo collapse are dictated by the 
potential, which, being a double integral over all space, is much rounder than the mass 
distribution. Therefore the effects of intrinsic triaxiality of initial density peaks are smaller 
than those due to the secondary perturbations, and so can be ignored. Furthermore, 
initial triaxiality is less severe in larger, 2-4 a peaks (Bardeen et al. 1986), which are the 
subject of the present study. With this caveat the smooth part of our density peak is still 
described by Eq. 

Moreover, in reality, the initial density peak will not be smooth, but will instead be 
sprinkled with many smaller scale positive and negative perturbations that arise from 
the same Gaussian random field that gave rise to the main peak. These secondary pertur- 
bations will perturb the motion of the dark matter particles from their otherwise purely 
radial orbits. 

So, in order to investigate effects such as tidal torques and nonradial collapse, it is 
necessary to consider the nonspherical portion of the density distribution. In addition 
to the smooth halo, RG87 and Williams et al. (2004) considered contributions from the 
secondary perturbations which arise from the same Gaussian random field that gave rise 
to the main halo. The overall initial density profile, linearly evolved to the present day, 
can be written as, 

p(x) = po[l + a (aO][l + eo(x)], (13) 

where po is the present day background density density excess due to the main halo 
is Sq(x), and is assumed to be spherically symmetric, and eo(x) is the density excess 
contributed by the random secondary perturbations. In a statistical sense, the growth 
rate will depend on So/ So, or, equivalently, x, either of which can be used to parameterize 
the strength of the tidal field in a spherically symmetric halo. The final expression for 
the growth rate is given by (Eq. [28] of RG87 and Williams et al. 2004) 

e ( x 0) = f^M (14) 

where eo (x) is the amplitude of the initial perturbation, given by 

eo(x) = ^Jd 3 k e k e* kx , (|e k | 2 ) = P(k) (15) 
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and fi and f 2 are functions of the 'time' parameter 8: fi(0) = 16 — 16 cos + sin 2 — 
99 sin 9, and f 2 (0) = 12 - 12 cos 9 + 3 sin 2 6 - 99 sin 9. 

At any given time and proper position the random acceleration due to perturbation 
field e(r) is given by (RG87 Eq. [40]): 

g(r,t)= gtot (r,t)-g 6 (r,t)«G / d 3 r , Pb{ f; t)e ^ t) (v'-r), (16) 



which is an integral over all space, and pj,(r,t) = po(t)[l + S(r,t)] is the background 
density due to the main halo, and is related to the total density given by Eq. [13] by, 
p(r,t) — pb(r,t)[l + e(r, t)]; po(t) is the average density of the Universe at epoch t. 
Remembering that the density distribution in the main peak and perturbation field and 
growth rate of perturbations are functions of radial position only, we use x instead of x, 
and r instead of r. Since, our goal is the rms value of acceleration, so scalar g will replace 
vector g. A major simplification is accomplished by decoupling the time dependence of 
acceleration, i.e. the rate of growth of acceleration, from its spatial variation. Let the 
initial acceleration field due to secondary perturbations be denoted by go(x). Then the 
dimcnsionlcss rate of growth of acceleration is given by, 

F g (x,t) = g{x,t)/g (x) = g{r,t)/g (r). (17) 
With these, the proper displacement of a particle can be evaluated as 



d p (x,t)= / dti / dt 2 g(x,t 2 ) = g (x) / dt x I dt 2 F g (x,t 2 ) = Ag (x)t F r (x,t), (18) 

JO Jo Jo JO 

Next we describe the two functions, F g (x,t) and go(x), separately. 



The acceleration growth rate becomes (RG87 Eq. [44]3) 

Fg (x,0) = ^ = *[i + '(M)MMMM) = 8 - m 

9o(x) poeox [/l( 0)_gg2 /2( 0)p 

The spatial dependence of acceleration is given by: 

Ag (d,x)=4Gp [[ P s (k)e- kd (l-^^)dk] 1/2 . (20) 

The comoving displacement is given by: 

d(x,t) = dp(x,t)/a(t) = Ag (x)tlF r (x,t)/a(t). (21) 

In the early part of the evolution of the halo most of the shells are still expanding. 
During this time secondary perturbations grow, and so does the acceleration, the velocity, 
and the displacement contributed by these perturbations to the particles in the shell. 
Because the secondary peaks are randomly distributed within the halo, they displace 
the dark matter particles in random directions from their original positions. This can be 
visualized as a shell having an internal velocity dispersion, resulting in a "swollen" shell. 
If we concentrate on a single particle, its orbit, viewed from the rest frame of the parent 
shell, will oscillate between an inner and an outer radius of that shell (i.e. pericenter, r p 
and apocenter, r a ). In general, this orbit will not be closed and will resemble a rosette. 
In a real situation, the "swell" of any given shell will gradually increase as the shell 
expands away from the center of the halo, and the influence of the secondary peaks 
grows, but for simplicity, as long as a given shell is expanding its dark matter particle 
positions and velocities are not corrected for the effects of secondary perturbations. So, 
the contributions are evaluated analytically but are not imparted to the shell until it 
reaches turnaround. 

The magnitude of the extra velocity imparted to a typical dark matter particle at the 
time of turn-around (Eq. [48] of RG87), 

|Av rms (M c /2)| = F v (x,t c /2)Ag [d(x,t c /2),x]t , (22) 



Note a typo in their paper: fi(0) in the numerator of Eq. (44) should be /2(f). 
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where the velocity growth factor F v (x, t) is: 
f 9 dt 

F v (x,t) = / F g {x,9 1 )—d6 1 (23) 

and the tangential, Vt an , and radial, v ra d, components of velocity are: 

2 1 
{Avtanf = -|Av rms (x)| 2 , and (Av rad ) 2 = -|Av rms (a;)| 2 . (24) 

The mean radial velocity at t — t c /2 is zero, while the mean tangential velocity is given 
by: 

2X l/2 



Mx) = [o) r m Av(x,t c /2) (25) 

It is important to note that j$ is inversely correlated with the height of the peak in 5, 
since both Aw and r m (the maximum expansion radius) decrease with an increase in 
the height v of the peak. The random angular momentum is important for the particle 
(shell) orbits since the larger it is, the larger is the orbital ellipticity, and then the shell 
penetrates less to the center, resulting in a flattening of the inner density profile. In 
several of the previously quoted papers (even Williams et al. 2004) only the random 
angular momentum was taken into consideration, while we take also account of ordered 
angular momentum, dynamical friction and adiabatic contraction. 

The collapse starts from the inner most shell, the one adjacent to the core. When 
the first shell reaches its r m it collapses and finds its r a and r p within the overall halo 
potential. It is assumed that the potential is changing slowly compared to the dynamical 
timescales of the shells, so that every shell conserves its adiabatic invariants, the radial 
and tangential momenta, 

jg(x) = Av tan r m (26) 

jr(x) = / V rad dr (27) 



throughout the collapse. This is an important assumption in the RG87 formalism, it is 
crucial to the computation of dynamics of shell crossing. Up to the moment when a given 
shell reaches its maximum expansion radius r m (x) at a time t = t c {x)/2 corresponding 
to 9 = 7T, the shell is assumed to be thin, its radial extent determined by the initial shell 
separation. At r m , the average dark matter particle in the shell is given its additional 
random velocity, Eq. [2U and Eq. [M] In either case the apocenter and pericenter can be 
calculated from the particle's energy integral and Aw ttm and Av rac i- The energy integral, 



E = iP(r m ) - \ 



, x 2 / . \ 2 
«M , I 38 

dt 



(28) 



is conserved if ip is not too rapidly changing with time, an approximation which is good 
outside the immediate region of the core. 
The radial velocity of a particle is 

v rad =[2(E-i;(r))-(Je/r) 2 ] 1 / 2 , (29) 

(RG87). Note that RG87 define their energy integral and potential as the negatives of 
the conventional definitions of these quantities. We use the conventional definitions, i.e. 
potential is a negative quantity inside the halo, and energy integral of a bound particle 
is negative. 

The addition of a random radial velocity ensures that the apocenter distance is greater 
than r m , the radius of maximum expansion in the absence of velocity perturbations. A 
measure of the mean radial momentum of the particle is the radial action, (RG87): 

j r = / v r dr = v r dr(r a — r p ) (30) 
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The radial distribution of mass within a shell between apocenter and pericenter is 
not uniform. The density in the radial range dr around r is proportional to the amount 
of time the particle spends there, (Eq. [53] of RG87), 

P(r)dr= r ! r "lf . ■ (31) 

ir p V rad dr 

If M s heii is the mass of the added shell, then the total mass distribution of the core 
and shell together is: 

Mi(r) - M(r) + M shM ( P(r')dr' (32) 



From this mass function, the new potential ip\{r) is calculated. The next mass shell is 
added and its probability distribution is calculated in the potential V , i( r )- However, the 
mass distribution of a newly added shell overlaps with shells which have fallen earlier; 
that is, the pericenter of the shell is at a smaller radius than the apocenters of some 
fraction of the previously added shells. Thus we must recompute after adding each new 
shell the orbits for each shell with which it overlaps. Since the potential is not violently 
changing, the adiabatic invariants of the orbits are conserved. In this case, the adiabatic 
invariants are the angular momentum jg and the radial action j r . By repeatedly adding 
shells in this manner, while adjusting the orbits so that the angular momentum and 
radial action of the orbits are conserved, a self-consistent mass distribution is built up. 

The collapse of a perturbation taking into account angular momentum, and dynamical 
friction may be calculated by solving the equation for the radial acceleration (Kashlinsky 
1986, 1987; Antonuccio-Delogu & Colafrancesco 1994; Peebles 1993): 

dv r h 2 (r, v)+j 2 (r, v) dr 

-dt = J5 ° {r) ~dt (33) 

where h{r, v) is the ordered specific angular momentum generated by tidal torques, 
j{r,v) the random angular momentum (see RG87), G(r) the acceleration, and /x the 
coefficient of dynamical friction. 

In the peculiar case of /i = 0, Eq. (f33|) can be integrated to obtain the square of 
velocity: 



v(rY = 2 



e - G I —^ d y + / ^ d v 

o y Jo y . 



m T (y) , , f r h 2 



(34) 



where e is the specific binding energy of the shell that can be obtained from the previous 
equation at turn-around when, dr / dt = 0. 

If /j, ^ 0, the previous equation must be substituted with: 



dv 2 „ 9 
— + 2^v 2 = 2 
dt 



h 2 + j 2 _ 



(35) 



which can be solved numerically for v. In the previous equation, we have two unknown 
quantities, the specific ordered angular momentum, h, and the coefficient of dynamical 
friction, fi. 

The ordered angular momentum is connected to the tidal interaction of the proto- 
structure with the neighboring ones. The explanation of galaxies spins gain through tidal 
torques was pioneered by Hoyle (1949). Peebles (1969) performed the first detailed cal- 
culation of the acquisition of angular momentum in the early stages of protogalactic evo- 
lution. More recent analytic computations (White 1984, Hoffman 1986, R88; Eisenstein 
& Loeb 1995; Catelan & Theuns 1996) and numerical simulations (Barnes & Efstathiou 



4 In other words, the mass particles, once position and velocities are assigned, are allowed to 
follow the appropriate orbit in the gravitational potential of the previously collapsed matter. As 
each mass shell is added, those previously added shells with which it overlaps have their orbits 
adjusted so that the angular momentum jg and radial moments j r integrated from pericenter 
to apocenter, are conserved. 
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1987) have re-investigated the role of tidal torques in originating the angular momenta 
of galaxies. 

In the present paper, we take into account both types of angular momentum: random 
j, and ordered, h. As described in RG87, to calculate the ordered angular momentum, 
one has first to obtain the rms torque, r(r), on a mass shell and then calculate the total 
specific angular momentum, h(r, v), acquired during expansion by integrating the torque 
over time (Ryden 1988a (hereafter R88), Eq. 35): 

V ' 7 3 W M sh ° J (^-sin< /3 /iW~/2W| V ' 

where M s h — 47rpb [1 + S(x)] x 2 Sx is the mass in a thin spherical shell of internal radius 
x, t is the tidal torque at time to, while the mean over-density inside the shell, S(r), is 
given by Eq. (fTTj) . 

As remarked by Del Popolo et al. (2001) the angular momentum obtained from Eq. 
(l36|) is evaluated at the time of maximum expansion i m . With the BBKS power spectrum 
(filtered on a galactic scale), for a v = 3 peak of mass ~ 2 x 10 11 M Q , the model gives a 
value of h = 2.5 x 10 74 

s 

It is interesting to compare this result with a different method like that of Catelan 
& Theuns (1996), who calculated the angular momentum at maximum expansion time 
(see their Eqs. (31)-(32)) and compared it with previous theoretical and observational 
estimates. Assuming the same value of mass and v used to obtain our previously quoted 
result and the same distribution of angular momentum as adopted by Catelan & Theuns 
(1996), we obtain a value for the angular momentum of h = 2.0 x 10 74 gc " in agreement 
with our result. 

The total angular momentum of a system is often expressed in terms of the dimcn- 
sionless spin parameter, 

A - B^l (37) 

where L is the angular momentum, summed over shells, and E is the binding energy of 
the halo. Expressing the values of angular momentum calculated as previously reported in 
terms of the spin parameter, we get a value of A ~ 0.04 for a 3a peak of ~ 10 12 AfeH. This 
is comparable with values found in the literature (Barnes & Efstathiou 1987; Vivtska et 
al. 2002). Using the parameters in Vivitska et al. (2002), the maximum of the distribution 
of A (well approximated by a log-normal distribution) is A = 0.035 while there is a 90% 
probability that A is in the range 0.02-0.1. 

Since our paper deals also with LSB galaxies, we must recall that LSB galaxies are 
more angular momentum dominated compared to normal galaxies with the same lumi- 
nosity (McGaugh & de Blok 1998). These galaxies are characterized by high values of A: 
as shown in Boissier et al. (2002), 35% of all galaxies (in number) having 0.06 < A < 0.21 
are LSBs. Therefore one expects according to several analytical papers (Nusser 2001; 
Hiotelis 2002; Le Delliou & Henriksen 2003; Ascasibar et al. 2004; Williams et al. 2004) 
that these objects should typically have shallower density cusps. Moreover the cores of 
these galaxies are much less dense than what the simulations indicate. In the case of 
LSB galaxies, we assumed a value of A = 0.06 performing a conservative estimate. In 
agreement with what is reported, larger values of A should imply a further flattening of 
the rotation curves in Fig. 3. 

A more complete description of this calculation can be found in Del Popolo (2006). 

It is important to note that after turn-around, not only the random angular momen- 
tum, j, will contribute to originate non-radial motions in the protostructure but also the 
ordered angular momentum, h. In fact, as successive shells expand to their maximum 
radius they acquire angular momentum and then fall in on orbits determined by the 
angular momentum. Actual peaks are not spherical; thus the infall of matter will not be 
purely radial. Random substructure in the region surrounding the peak will divert in- 
falling matter onto non-radial orbits. The role of these random motions is of fundamental 
importance in structure formation. 

5 Notice that there is a very mild dependence of A on peak height, or mass 
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We took into account dynamical friction by introducing the dynamical friction force 
in the equation of motions. In ordered to calculate fj,, we recall that in hierarchical 
universes, matter is concentrated in lumps, and the lumps into groups and so on, which 
act as gravitational field generators. One can calculate the stochastic force generated by 
these field generators and then following Kandrup (1980) the dynamical friction force 
per unit mass 

- ^te^ log (I.!***.) -iL . _ A _iL (38) 

where N = ^-R^ ys n a , is the total number of field particles, R sys the system radius, m a 
and n a are respectively the average mass and the number density of the field particles, 
n ac = n a x a 3 is the comoving number of field particles and a the expansion parameter, 
connected to the proper radius of a shell through: 

r(n,t) =na(n,t) (39) 

The number and mass of the field generators is then calculated using the theory 
of Gaussian random fields (BBKS). The acceleration obtained from Eq. (f3"5|) is used in 
the equation of motion. Again, the reader is referred to Del Popolo (2006) for a more 
complete description. 

The shape of the central density profile is influenced by baryonic collapse: baryons 
drag dark matter in, causing the so called adiabatic contraction (AC) steepening the 
dark matter density slope. By considering a spherically symmetric protostructure, before 
dissipation, consisting of a baryonic fraction F with F << 1 of dissipational baryons and 
a fraction 1 — F of dissipationless dark matter particles constituting the halo, Blumenthal 
et al. (1986) described an approximate analytical model to calculate the effects of AC. One 
assumes that the dissipational baryons and the halo particles are well mixed initially (i.e., 
the ratio of their densities is F through the protostructure). This is due to the fact that 
the original angular momentum of the dark matter halo comes from gravitational (tidal) 
interactions with its environment. Thus, the dark matter and the gas experience the 
same torque in the process of halo assembly and should initially have (almost) the same 
specific angular momentum (Klypin et al. 2002). So, a usual assumption is that initially 
baryons had the same density profile as the dark matter (see the previous discussion and 
Mo et al. 1998; Cardone & Sereno 2005; Treu & Koopmans 2002; Keeton 2001). 

In real systems, baryons assume a distribution in the central regions of the halo 
which is determined by the competition between dissipation and star formation. The 
dissipative infall of baryons will cease either when most of the baryon gas is converted 
to stars, ending dissipation, or when a rotationally supported gas disk is formed. As 
the baryons dissipate energy, they fall toward the center of the halo. If the baryons 
conserve their angular momentum as they fall inward, then the baryons will dissipate 
energy until they form a rotationally supported disk, whose shape is determined by the 
angular momentum distribution h(M) of the baryons. So, summarizing, the two methods 
of ending dissipative infall result in a spheroidal distribution (if star formation ends the 
infall) or in a disk of stars (if angular momentum ends the infall) (R88a) . So following the 
usual practice, the final baryon distribution is assumed to be a disk (for spiral galaxies) 
(Blumenthal et al. 1986; Flores et al. 1993; Mo et al. 1998; Klypin et al. 2002; Cardone 
& Sereno 2005). In our calculations, we shall assume the Klypin et al (2002) model for 
the baryon distribution (see their subsection 2.1), when dealing with mass scales typical 
of spiral galaxies. In the case of elliptical galaxies and clusters a typical assumption (that 
we shall use in the paper) is that baryons collapse to a Hernquist configuration (Rix et 
al. 1997; Keeton 2001; Treu & Koopmans 2002). One objection that one could advance 
is that clearly if the final baryon distribution is assumed to be a disk, the inner gravity 
potential is clearly non-spherical. The usual assumption that is made at this step (see 
RG87; R88a,b; Mo et al. 1998; Klypin et al. 2001; Cardone et al. 2006) is that the final 
configuration (in this case the disk) is assembled slowly, so that we can assume that 
the halo responds adiabatically to the modification of the gravitational potential and it 
remains spherical while contracting. The angular momentum of the dark matter particles 
is then conserved and a particle which is initially at radius r ends up at a rdius r', as 
described in the following. 
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The scale of the Hernquist profile is fixed by the competition between dissipation and 
star formation, when most of the baryon gas is converted to stars (see RG87; Gncdin et 
al. 2004). 

We take into account adiabatic compression as done in RG87 improving it as done 
by Gnedin et al. (2004) (sec the following). The Blumcnthal ct al. (1986), RG87, R88ab 
model for the adiabatic contraction can be described as follows. Since the baryon fraction 
F is much less than one, the adiabatic invariants of the dark matter orbits are conserved 
(Blumcnthal et al. 1986; RG87). If the matter is on circular orbits, the invariant is 
rM(r); if the dark matter is on radial orbits in a self-similar distribution, the invariant is 
r m M(r m ) (Blumcnthal et al. 1986). Since the actual orbits are neither circular non radial, 
the two adiabatic invariants which are preserved are the angular momentum jg and j r , 
of the orbits (RG87). In the absence of dissipation, the mass distribution Mo(r) has a 
value of F which is constant throughout the structure. A self-consistent final solution, 
in which the baryons form a disk (or Hernquist configuration), conserving h, and the 
dark matter is compressed, conserving jg and j r , is found by an iterative procedure. 
In the zeroth-order approximation, the baryons have their predissipation distribution 
M B o = FM {r), and the dark matter has the distribution M D0 = (1 - F)M (r). While 
the dark matter distribution is held constant, the baryons fall inward, preserving their 
angular momentum until they are on circular orbits. Baryons originally at radius r and 
with specific angular momentum h will end up, in the stiff halo, at a radius, 



The baryons now have a new, more compact distribution Mgi(r), and the central con- 
centration of the baryons will draw the dark matter inward. The new potential of the 
baryons plus the dark matter is given by, 



After fixing the value of jg for each mass shell of the halo, the value of the apocenter is 
adjusted until the orbits of the dark matter in the potential <&i(r) have the same value of 
j r , which they had in the predissipation distribution. The distribution of mass Mm(r) 
of the dark matter in the potential <I>i is built orbit by orbit, ensuring that jg and j r are 
conserved. In the second iteration one has: 



and the new potential $2 is calculated and the dark matter orbits are adjusted to preserve 
]g and j r . The iteration continues till convergence. The previous model can be improved in 
two ways, as shown by Gnedin et al. (2004) and Klypin et al. (2002). It is worth stressing 
that there is some debate of the validity about the adiabatic compression formalism, with 
authors like Jesseit et al. (2002) finding a substantial agreement between the final dark 
matter distribution in numerically simulated haloes and that predicted by the adiabatic 
compression approach. On the other hand, this result has been contradicted on the basis 
of a set of higher resolution numerical simulations by Gnedin et al. (2004). According to 
these authors, the standard adiabatic compression formalism systematically overpredicts 
the dark matter density profile in the inner 5% of the virial radius, and the adiabatic 
compression formalism overpredicts the dark matter density less than ~ 10% at r/r v 
~ 0.1, while the error quickly decreases for larger values of r/r v . Gustafsson et al. (2006) 
confirmed that the Blumcnthal ct al. (1986) model overestimates the central dark matter 
density. Moreover they showed that the modified model proposed by Gnedin et al. (2004) 
even if it is a considerable improvement it is not perfect. Moreover it is found that 
the contraction parameters in their model not only depend on the orbital structure of 
the dark-matter-only halos, but also on the stellar feedback prescription which is most 
relevant for the baryonic distribution. As a caveat we want to recall that according to 
Romano-Diaz et al. (2008), Gustafsson et al. (2006), even if focused on the AC, had 
an insufficient resolution to reach conclusive results. In the following, we shall take into 



G[M B0 (r) + M D0 (r')Y 



(40) 




(41) 



G[M Bl (r) + M D1 (r")} 



(42) 
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account the Gnedin et al. (2004) simple modification which describes numerical results 
more accurately than Blumenthal et al. (2006). They proposed a modified adiabatic 
contraction model based on conservation of the product of the current radius and the 
mass enclosed within the orbit-averaged radius, 



where T r is the radial period, r a is the apocenter radius and r p the pericenter radius. 
The previous, classical AC model, assumes no angular momentum exchange between 
different components (e.g., baryons and dark matter), transfer of angular momentum 
from the baryons to the dark matter can be produced by dynamical friction, and the 
change of the AC model due to this effect can be taken into account by means of a 
simple model by Klypin et al. (2002). 

For what concerns the universal baryon fraction, it has been inferred by observa- 
tions involving different physical processes (Turner 2002). The power spectrum of matter 
inhomogeneities of large-scale structure is sensitive to the quoted F; the Two Degree 
Field Galaxy Redshift Survey has reported a value of 0.15 ± 0.07 (Percival et al. 2001). 
Measurements of the angular power spectrum of the CMBR also provide a very sig- 
nificant estimate. The combined analysis in Jaffe et al. (2001) of several data gives 
F = 0.186tQ'QQg. By using Wmap collaboration data (Spergel et al. 2003) one obtains 
F ~ 0.16 = 0.047 ± 0.006; n m = 0.29 ± 0.07). After baryons cools and form stars, 
the baryon-to-total mass at the virial radius does not need to be equal to the previ- 
ous quoted universal baryon fraction and may deviate from it depending on feedback 
effects, hierarchical formation details, and heating by the extragalactic UVB flux. In the 
AC calculation one has to take account that for large galaxies, only about half of the 
baryons, in principle available within the virial radius, goes into the central galaxy. For 
dwarf galaxies, this amount may be much smaller due to feedback effects, hierarchical 
formation details, and heating by the extragalactic UVB flux. In the following, in the 
case of large galaxies, we use one half of the universal baryonic fraction and for dwarfs 
values obtainable from Hoeft et al. (2007) (their Fig. 1). For example in the case of a 
dwarf galaxy of 1O 1O M it is about 0.06. 

3. Results and discussion 

After fixing the initial conditions and describing how to calculate angular momentum, 
dynamical friction and adiabatic contraction, we can use the model in Section 2 to obtain 
the density profile of haloes. In Fig. 1, we plot the profiles obtained through our model and 
those predicted by numerical simulations. In Fig. 1, the short-dashed line, the dashed line, 
and the solid line represents, respectively, the density profile for haloes of lO 12 /!" 1 M©, 
10 10 /i _1 M©, and 10 8 h~ 1 M Q calculated by means of the model of this paper. In order to 
compare the results of the model to those of N-body simulations, we plotted the NFW 
profile (dot-dashed line) for haloes having masses equal to lO 12 ^ 1 M©, and c = 10, 
and the results for the same mass obtained in the Aquarius Project by Navarro et al. 
(2008) (upper dashed- lines). Navarro et al. (2008) showed that density profiles deviate 
slightly but systematically from the NFW model, and are better approximated by a 
fitting formula where the logarithmic slope is a power-law of radius: the Einasto profile, 



where r_2 and p_2 are connected to the scaling radius and density of the NFW profile 
(see the following) by r_ 2 = r s and = ps/4. 

Density profiles become monotonically shallower inwards, down to the innermost re- 
solved point, with no indication that they approach power-law behavior. The innermost 
slope they measure is slightly shallower than -1, a result supported by estimates of the 
maximum possible asymptotic inner slope. Shallower cusps, such as the r~ - 75 behavior 



M(r)r — const, 



(43) 



where the orbit-averaged radius is 




(44) 



Hp{r)/p. 2 ) = -2/a[(r/r_ 2 ) Q - 1], 



(45) 
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predicted by the model of Taylor & Navarro (2001), cannot yet be excluded. The NFW 
profile for the given mass was calculated by means of the relationships connecting the 
concentration parameter, c, and the virial mass, M v , to the shape of NFW profile. We 
used the following equation for c, 



c - 13 - 6 U^J ' (46) 

(Gentile et al. 2007), and the usual one for the NFW profile, 

p(r) = P -± = ^ (47) 

P{ ' r/r s {l + r/r s f r/r s (l + r/r s )* 1 ' 

where 

A„ c 3 

= "Tio g (i + c )- c /(i + c )' (48) 

and A„ is the virial overdensity (see Bryan & Norman 1998). The scaling radius, r s , of 

the NFW profile is connected to the virial radius, concentration parameter and virial 
mass through c — r v /r s , where, 



10 n M P 



kpc (49) 



(Gentile et al. 2007). The residuals from the best Einasto fits are extremely small, and 
show no sign of systematic deviations down to the innermost resolved point. NFW resid- 
uals are then less than 30 % over the full radial range and below 20% within r_ 2 . 

Figure 1 shows that halos generated using our method are different in character 
from the profiles predicted by numerical simulations, like those of NFW or the Aquarius 
Project. This result is not surprising since the types of evolution that numerical N-body 
and our halos undergo are rather different: the former are produced as a cumulative result 
of many minor and major mergers of smaller sub-halos, while the latter are the product 
of quiescent accretion of lumpy matter onto the primary halo. The NFW profiles change 
slope rapidly from a = 1 to a = 3 at the characteristic radius r s . For example, in the case 
of the lO 12 ft. -1 -M0, the NFW profile have a characteristic scale length, equal to 0.1r„, 
beyond which the density profile steepens, so that much of the mass is piled up within 10% 
of the virial radius. It is interesting to note that Navarro et al. (2004), Stadel et al. (2008) 
and the Aquarius Experiment start to show that the spherically-averaged density profile 
is not described, as in the NFW model just by an outer power-law p oc r~ 3 and inner 
power-law p a r 1 but becomes progressively shallower inwards, and at the innermost 
resolved radius, the logarithmic slope is 7 = — din p/dlnr ^ 1, ruling out claims of a 
steep p oc r -1,2 central cusp, and not excluding shallower cusps as the r -0 75 behavior 
predicted by the model of Taylor and Navarro (2001). In other words, simulations with 
higher resolution are merging to the conclusion that the inner slope is not so steep as 
it appeared in less recent simulations and so in some cases even marginally consistent 
the DM halo of NGC 4605 (p oc r" 65 ) (Simon et al. 2003). In the following, we discuss 
another reason why density profiles obtained in N-body simulations are different from 
those obtained in our model. 

NFW profiles change slope rapidly from a = 1 to a = 3 at the characteristic radius r s . 
For example, in the case of the 10 12 /i _1 M Q model, the NFW profile has a characteristic 
scale length, equal to 0.1r„, beyond which the density profile steepens, so that much 
of the mass is piled up within 10% of the virial radius. The haloes obtained using the 
model in Section 2 are different in character from the profiles predicted by numerical 
simulations, like those of NFW. Within the virial radius the log-log density slope changes 
gradually and the slopes of the inner part of haloes flatten with decreasing mass. Our 
results show a steepening of the density profile with increasing mass with slopes a ~ 
for M ~ 10 s — 10 w h- 1 M Q , and a ~ 0.8 for M ~ 1O 12 /i _1 M . The profiles with M ~ 
10 8 — 10 10 /i _1 M Q are well fitted by means of a Burkert's profile considered a good fit 
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to the dark matter rotation curves inferred from observations (e.g., Salucci & Burkert 
2000). The functional form of this profile is characterized by, 

P(r) = (l + r/r )[i°+(r/r )T (50) 

where p Q ~ p s and r Q ~ r s (EZ01). The dark matter within the core is given by M = 
1.6p rg. Although the dark matter parameters r a , p and M Q are in principle independent, 
the observations reveal a clear connection (Burkert 1995): 

M a = 4.3 x 10 7 {r /kpc) 7/3 M Q . (51) 

In panel (a) of Fig. 1, we plot the result of the model of this paper in the case 
baryonic infall is not taken into account, while in panel (b) baryonic infall is taken into 
account. Comparing panel (a) and (b), we observe that the presence of baryons leads to 
a steeper profile in agreement with previous studies (Blumenthal et al. 1986; Flores et 
al. 1993; Klypin et al. 2002; Gnedin et al. 2004; Gustafsson et al. 2006). It is important 
to note that the effect of baryonic infall is more evident at early redshifts, as shown in 
Fig. 5. As shown in Fig. 5, from z = 3 to z = 2 the profile becomes even steeper than 
the NFW cusp. Later the effect of angular momentum and dynamical friction reduces 
the slope and erases the cusp. However, as shown in Fig. 1 the initial steepening of the 
profile remains visible in the profiles calculated through the model of the present paper 
(M ~ HP/i^Mq, M ~ 10 w h- 1 M Q , and M ~ lO^^M©), since the profile which 
takes into account baryonic infall (panel b) are steeper than those not taking it into 
account (panel a). 



The model of the present paper gives just one profile for a given halo mass. In reality 
N-body simulations and analytical models different from those of Ryden- Williams and 
that of this paper (e.g. the extended Press-Schechter formalism), predict an ensemble of 
profiles for a given halo mass, depending on the different formation histories. For example, 
using the extended Press-Schechter (EPS) conditional probabilities for halo progenitor 
masses one can construct detailed histories of the mass assembly of dark matter haloes. 
One is interested in computing the mass accretion histories (hereafter MAHs) of dark 
matter haloes, defined as ^(M ,z) = M(z)/M , where M(z) is defined as the mass of 
the main progenitor halo, and M is the halo mass at z = 0. One then defines the average 
mass accretion history of a halo of mass M as, 

N 

*(M ,z) = l/Nj2*i(M ,z), (52) 
»=i 

where the summation is over an ensemble of N random realizations (see van den Bosch 
& Frank 2002). One then can convert the mass growth curves M(z) to density profiles 
as for example shown by Nusser & Sheth (1999). The single profile that we obtain in the 
present paper must be considered as the average profile of simulations. 

Two important things must be noticed: a) less massive haloes are less concentrated; 
b) the halo's inner slope is smaller for smaller mass. 



The first point can be explained as follows: higher peaks (larger u), which are progen- 
itors of more massive haloes, have greater density contrast at their center, and so shells 
do not expand far before beginning to collapse. This reduces j and h and allows haloes 
to become more concentrated. An alternative explanation is connected to the quoted 
angular momentum-density anti-correlation showed by Hoffman (1986): j oc jy~ 3 / 2 (and 
similarly for h). So, density peaks having low (high) value of v acquire a larger (smaller) 
angular momentum than high v peaks and consequently the halo will be less (more) 
concentrated. It is important to notice that the quoted trend of increased central concen- 
tration as a function of mass applies only to halos that started out as peaks in the density 
field smoothed with a fixed Rf scale. Our conclusions do not mean that, for example, 
clusters of galaxies will be very much more centrally concentrated than galaxies, since 
different smoothing scales would apply in the two cases. 



16 



A. Del Popolo: Density profiles of dark matter haloes 



Point (b) can be explained in a similar way to (a), as described previously. Less mas- 
sive objects are generated by peaks with smaller which acquire more angular momen- 
tum (h and j). The angular momentum sets the shape of the density profile at the inner 
regions. For pure radial orbits, the core is dominated by particles from the outer shells. As 
the angular momentum increases, these particles remains closer to the maximum radius, 
resulting in a shallower density profile. Particles with smaller angular momentum will be 
able to enter the core but with a reduced radial velocity compared with the purely radial 
SIM. For some particles the angular momentum is so large that they will never fall into 
the core (their rotational kinetic energy makes them unbound). Summarizing, particles 
with larger angular momenta are prevented from coming close to the halo's center and 
so contributing to the central density. This has the effect of flattening the density profile. 
This result is in agreement with the previrialization conjecture (Peebles & Groth 1976; 
Davis & Peebles 1977; Peebles 1990), according to which initial asphericities and tidal 
interactions between neighboring density fluctuations induce significant non-radial mo- 
tions that oppose the collapse. In order to reproduce the NFW profile, we performed an 
experiment similar to that performed by Williams et al. (2004), namely we reduced the 
magnitude of the h and j angular momentum, dynamical friction, [i and we considered 
the system as constituted only of dark matter {F = Q). The experiment was performed on 
the halo of mass 1O 12 /i -1 M0, and in order to reproduce the NFW profile having c = 10 
and mass ~ 10 12 h~ 1 M Q , we had to reduce the magnitude of h by a factor of 2, j and [i 
by a factor of 2.5. The result of the quoted experiment is the dashed line in Fig. 2, which 
closely reproduces the NFW profile (dot-dashed line), and the dotted line is the profile 
for the 10 12 h~ 1 MQ halo obtained by means of our modelH. 

Similarly, Williams et al. (2004) had to reduce random velocities, which amount to 
reducing the angular momentum, in order to obtain a NFW profile. With each reduction 
of the random velocities, the profiles get steeper at the center. This effect can be un- 
derstood, as already reported, as follows: the central density is built up by shells whose 
pericenters are very close to the center of the halo. Particles with larger angular momenta 
are prevented from coming close to the halo's center and so contributing to the central 
density. The correlation between increasing angular momentum and the reduction of in- 
ner slopes in halos has been also noticed by several other authors (Avila-Reese et al. 1998, 
2001; Subramanian et al. 1999; Nusser 2001; Hiotelis 2002; Le Delliou & Henriksen 2003; 
Ascasibar et al. 2004). As previously mentioned, there could be, according to Williams 
et al. (2004) another reason of the difference of density profiles obtained from N-body 
simulations and models like that of this paper (or William's). As shown in their Fig. 1 
(lower-left panel), the specific angular momentum distribution (SAM) in NFW-like halo 
is more centrally concentrated than the SAM distribution of the reference halo obtained 
with their model, model which is similar to that of this paper, and is closer to those 
of typical halos emerging from numerical simulations. This may suggest that haloes in 
N-body simulations lose a considerable amount of angular momentum between 0.1 and 1 
r v . Since virialization proceeds from inside out, this means that the angular momentum 
loss takes place during the later stages of the halos' evolution, rather then very early on. 
This is somehow confirmed by the so called angular momentum catastrophe, namely the 
fact that dark matter halos generated through gas-dynamical simulations are too small 
and have too little angular momentum compared to the halos of real disk galaxies, pos- 
sibly because it was lost during repeated collisions through dynamical friction or other 
mechanisms (van den Bosch et al. 2002; Navarro & Steinmetz 2000). The problem can 
be solved invoking stellar feedback processes (Weil et al. 1998), but part of the angular 
momentum problem seems due to numerical effects, most likely related to the shock cap- 
turing, artificial viscosity used in smoothed particle hydrodynamics (SPH) simulations 
(Sommer-Larsen &Dolgov 2001). 

We discussed the effect of changing the magnitude of angular momentum but we 
did not speak of the effect of changing the magnitude of /i (dynamical friction). The 
effect of changing this last quantity is very similar to changing the magnitude of angular 
momentum: an increase in the term \x produces shallower profiles as larger values of 



Note that the density profile of the Aquarius Experiment can be reproduced in the same 
manner with a small change in the h and j parameters used for the NFW profile. 



A. Del Popolo: Density profiles of dark matter haloes 



17 



angular momentum does. This is expected from Del Popolo (2006) (Fig. 1), showing that 
dynamical friction influence the dynamics of collapse in a similar way to that of angular 
momentum slowing down the collapse of outer shells and so compelling the particles to 
remain closer to the maximum radius. 

At this point, we have to stress that the flattening of the density profile is caused 
by two main effects: angular momentum and dynamical friction. The other effect that 
we have considered, namely baryonic infall, counteracts the effect of the density profile 
flattening produced by angular momentum and dynamical friction. The role of baryonic 
infall is predominant at early times and produces a steepening of the profile, while later 
its effect is overwhelmed by the effect of angular momentum and dynamical friction. The 
result is similar to the one described by Williams et al. (2004) with the difference that in 
the present paper the effect of ordered angular momentum and dynamical friction adds 
to that of random angular momentum studied by Williams et al. (2004) with the result 
that the flattening of the density profiles is larger. 

In Fig. 3, we plot the rotation curves obtained by our model and we compare them 
to four LSB galaxies shown in Williams et al. (2004). The four LSBs in Fig. 3 were 
taken from the high-resolution data of de Blok & Bosma (2002). Of the 26 LSBs that are 
presented in that paper we picked the ones with high inclination angles, smooth rotation 
curves, and good agreement between HI, Ha, and optical data. In all the four cases, the 
data are compared with the rotation curve obtained using our model (solid line) and 
with rotation curves obtained from NFW profile (dotted lines), given by, 

= ^ f Ml + CT ) -«/(! + ex) l 1 ^ 



x[]n(l + c) - c/(l + c)] 

where x = r/r v and V v is the virial velocity 0- The concentration parameter c was chosen 
to be consistent with the NFW predictions, calculating it from the mass of the galaxy 
through Eq. (I46| . For all the four galaxies the rotation curves obtained with the model 
of the present paper (solid lines) are a considerably better fit than NFW (dotted lines). 
Fig. 3 shows that NFW haloes are higher than the rotation curves obtained using our 
model, in which more massive haloes tend to be more centrally concentrated and have 
flatter rotation curves. Less massive haloes are less concentrated, and have slowly rising 
rotation curves. In contrast NFW rotation curves rise very steeply and as a consequence 
NFW fits to dwarf galaxy rotation curves have too low concentration parameters (van 
den Bosch & Swaters 2001) compared to N-body predictions. NFW fails to reproduce the 
velocities and the shape of the observed rotation curves, predicting too high velocities 
in the central part of haloes, and even leaving c as free parameter, instead of using Eq. 
(T46|) . there is no appreciable improvement in the fit. Using Eq. ([46]) . one obtains very low 
values of c. The result is similar to that described by Gentile et al. (2004): data are much 
better described by core-like profiles, like the Burkert profile generating flatter rotation 
curves, 

1 /2 

V(r) = ( 2 ^° r ° 3 )) {in [(1 + r/r )V(l + (r/r ) 2 )] - arctan(r/r )} 1/2 , (54) 

Our rotation curves are very similar to those generated by the Burkert profile and the 
residuals and discrepant points for our rotation curves are close to that given in Gentile 
et al. (2004) for the Burkert 's fit to their data. 

Our halos appear to be a closer match to the halos of spiral and dwarf galaxies, than 
are N-body halos. This may indicate that the halos of real late- type disk galaxies undergo 
a formation scenario similar to the one depicted by our method, i.e. collapse proceeds 
through a quiescent accretion of lumpy material and minor mergers, rather than through 
a merger-driven formation process characteristic of fully hierarchical models. 

Before going on, we notice that of the four rotation curves, NGC100 seems to be 
almost equally well fitted by the two type of models. In some other cases, the difference 
between the NFW rotation curves and those of the present paper is even smaller than 



7 The value of the characteristic velocity V v of the halo is defined in the same way as the virial 
radius r v . 
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for NGC100 and so it would be worthwhile to discuss halo-to-halo scatter. During the 
hierarchical assembly of dark matter haloes, the inner regions of early virialized objects 
often survive accretion onto a larger system, thus giving rise to a population of subhaloes. 
Depending on their orbits and their masses, these subhaloes therefore either merge, are 
disrupted or survive to the present day. To fully describe, in a statistical sense, the 
non-linear distribution of mass in the Universe, it is essential that halo substructure is 
taken into account. These subhalos appear ubiquitous in high resolution cosmological 
simulations and provide the source of fluctuations. One can study the evolution of halos 
under the influence of a generation of subhalos, while real halos grow continuously by 
accretion and mergers. Fluctuations due to subhalos in parent halos are important for 
understanding the time evolution of dark matter density profiles and the halo-to-halo 
scatter of the inner cusp seen in recent ultrahigh resolution cosmological simulations. This 
scatter may be explained by subhalo accretion histories: when we allow for a population 
of subhalos of varying concentration and mass, the total inner profile of dark matter can 
either steepen or flatten. 

As noticed in the introduction, while on galactic scales a large number of studies 
predicts central cores and it seems that the cusp/core problem is a real problem not 
attributable to systematic errors in the data (de Blok, Bosma & McGaugh 2003), on 
cluster scales the situation is less clear with slopes ranging from a — 0.5 (e.g. Sand et al. 
2002; Sand et al. 2004) to a = 1.9 (Arabadjis et al. 2002). 

In order to study the problem on cluster scales, we have calculated the density pro- 
file evolution of dark matter and that of the total matter distribution for a halo of 

We used the method of Sect. 2 to calculate the density profile and then we repeated 
the experiment shown in Sect. 3 to reproduce the NFW profile for the quoted mass that 
we considered as the initial condition, d 

After we evolved this profile using the method of Sect. 2 and taking the redshift 
dependence in the model using the technique described in Del Popolo (2001) (the reader 
is referred to the quoted paper to have more insights). The goal was to study the evolution 
of a NFW profile similarly to El-Zant et al. (2001,2004), Tonini et al. (2006), Romano- 
Diaz et al. (2008). 

Before showing the results, we recall that in real haloes angular momentum is lost in 
the final stages of collapse through the transfer of angular momentum from the subclumps 
to external material through dynamical friction (Quinn & Zurek 1988; R88b; Klypin et al. 
2002). Once the baryons condense to form stars and galaxies, they experience a dynamical 
friction force from the less massive dark matter particles as they move through the halo. 
Energy and angular momentum is transferred to dark matter, increasing its random 
motion. Moreover, angular momentum acquired in the expansion phase gives rise to non- 
radial motions in the collapse phase. The spherical approximation that we use, which 
ignores the possibility of large substructure forming, neglects the transfer of angular 
momentum through the interactions of subclumps q In the present paper, the effect 
of dynamical friction was just taken into account as an additive term in the equation 
of motion, Eq. (|33[) . without changing the sphericall symmetry of the problem as done 
by Peebles (1993) and Del Popolo (2006). The effect of the dynamical friction term in 
the equation of motion is very similar to that of angular momentum: an increase in the 
term \i produces shallower profiles as larger values of angular momentum does. This is 
expected from Fig. 1 of Del Popolo (2006), showing that dynamical friction influences 
the dynamics of collapse in a similar way to that of angular momentum slowing down the 
collapse of outer shells and so compelling the particles to remain closer to the maximum 
radius. 

Fig. 4a plots the evolution of a density profile of 1O 14 /i _1 M0. The solid line repre- 
sents the initial density profile at z = 3, which slightly steepens due to baryon settling 

8 We could have directly started from a NFW profile for the given mass calculated by means 
of the NFW fitting formula. 

9 An approximate analytical approach to the problem of exchange of angular momentum 
between the dark matter particles and infalling baryons has been developed by Klypin et al. 
(2002). 
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p(r) - — , n ,: b :\,^ 3 - a (55) 



in virialized dark matter haloes (AC) (dot-short-dashed line) at z = 2. Moreover, angu- 
lar momentum acquired in the expansion phase gives rise to non-radial motions in the 
collapse phase. The effect of angular momenta and dynamical friction overcomes that of 
the AC and the profile starts to flatten (z — 1.5 dotted line; z — 1 short-dashed line; 
z = long-dashed line). The final dark matter profile (long-dashed line) is characterized 
by a log-log slope of a ~ 0.6 at 0.01r s . So the situation is similar to that of haloes on 
galactic scales but the slope is larger than for dwarf galaxies. In Fig. 4b, we plot the 
evolution of the total density profile at initial redshift, z — 3, (solid line) and at the final 
redshift z = (dotted line). The plot shows that the total density profile slightly changes 
with time and the cusp is not "erased" as in the case of the dark matter profile. This 
result also implies that the baryonic component becomes steeper than the original NFW 
profile. The behavior of total mass is in agreement with X-ray observations by Chandra 
and XMM (Buotc 2003, 2004; Lewis et al. 2003), weak lensing (Dahle et al. 2003) and 
strong lensing (Bartelmann 2003), which are consistent with a cusp having a = 1 or 
larger. The behavior of the dark matter halo is in agreement with analysis of Sand et al. 
(2002, 2004) who fitted the baryonic and dark matter profiles only in the very inner part 
of the cluster MS 2137-23 within ~ 50ft -1 kpc by means of a generalized NFW profile: 

(r/r s ) a (l + r/rs) 

obtaining a nearly flat core a = 0.35. The steepening of the baryonic component is 
consistent with what found by Brunzendorf & Mcusingcr (1999), who found that the 
projected galaxy distribution in the Perseus cluster diverges as r . 

The results previously reported have several implications on the effort to test predic- 
tions of the CDM model observationally. The test that received much attention in the 
last decade, as several times stressed, is the density distribution in the inner regions of 
galaxies and clusters. 

On galactic scales, as our results show, the infall of baryons at early times steepen the 
cusp due to AC (in agreement with previous studies: Blumenthal et al. 1986; Gncdin et al. 
2004; Gustafsson et al. 2006) but later the cusp is erased through dynamical friction and 
non-radial motions effects. This results are in agreement with other analyses which study 
separately the effects of dynamical friction (e.g., EZ01) and those of angular momentum 
(non-radial motions) (e.g., Nusser 2001; Hiotelis 2002; Ascasibar et al. 2004; Williams et 
al. 2004) and with the recent Sph simulations of Romano-Diaz et al. (2008). 

Going to larger scales the situation changes. The analysis of the density distribution 
for the bright galaxies is complicated by the uncertain contribution of stars to the total 
mass profile (Treu & Koopmans 2002; Mamon & Lokas 2004). Some analyses tend to 
favor inner slopes shallower than predicted by CDM (e.g., Gentile et al. 2004) but others 
deduce slopes of the inner profiles that are at least marginally consistent with predictions 
(Treu & Koopmans 2002, 204; Koopmans & Treu 2003; Jimenez et al. 2003). As previously 
reported, our results shows a steepening of the density profile with increasing mass with a 
density profile of haloes of mass > 10 12 ft- _1 M Q having slopes > 0.8. This is in agreement 
with recent N-body simulations having a logarithmic slope that decreases inward more 
gradually than the NFW profile (Hayashi et al. 2003; Navarro et al. 2004; Stadel et al. 
2008). In the case of Stadel et al. (2008) the logarithmic slope is 0.8 at 0.05% of r v . 

The density distribution in clusters of galaxies can, in principle, provide a cleaner test 
of the models because the effects of the baryons and gas on the dark matter distribution 
are expected to be smaller and simpler. However also in this case observations predict 
slopes ranging from a < 0.5 (a = 0.35, Sand et al. 2002, 2004) to values larger than one 
(Arabadjis et al. 2002) with in some cases different results even for the same object. This 
is the case of the cluster MS2137-23 studied by Sand et al. (2004), who found a shallow 
density slope (< 0.5) while Dalai & Kecton (2003), Bartelman & Meneghetti (2004) and 
Gavazzi (2003) contested Sand's results, which according to them is neglecting lens ellip- 
ticity, and found consistency of the inner slope with a NFW profile. Our result concerning 
cluster scales makes a difference between dark matter and total mass distribution: the 
first tends to be less cuspy than the NFW profile in agreement with some observations 
(e.g., Sand et al. 2002, 2004), while the second is a bit more cuspy than the NFW profile 
(in agreement with Brunzendorf & Meusinger 1999) . 
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In the present paper, we did not take into account the baryon-DM interaction. As 
previously reported once the baryons condense to form stars and galaxies, energy and 
angular momentum is transferred to dark matter, increasing its random motion. 

In future work, we will develop a semianalytic model to follow the evolution of the 
baryonic component, and its interaction with DM, and we will compare these with full 
hydrodynamical simulations. This should bring new insights into galaxy formation, and 
directly address possible small-scale challenges to the CDM theory. The role of the inter- 
action between DM and baryons could bring further changes in the final density profile, 
that we shortly discuss. It is interesting to note that the problems of CDM only become 
clear on length scales where the baryons start playing a role and that this applies not 
only to the cusp/core problem but also to the missing dwarfs problem. Moreover, it is 
noteworthy that if, indeed, most starforming galaxies in the early universe lost their DM 
cusps because of stellar feedback, the missing dwarfs (satellites) problem could also be 
solved. Dwarf galaxies without a central cusp have a lower average core density than 
cuspy ones, and are hence much easier to disrupt tidally during the hierarchical assembly 
of larger galaxies (Mashchenko & Sills 2005). As a consequence, the removal of galactic 
cusps by stellar feedback in the early universe would result in fewer satellites today. This 
again indicates that baryon physics could be one of the missing pieces of the puzzle, and 
will very likely make a major contribution toward a solution. If this is true, it would be 
unwise to ignore the conclusions to which data arc leading us, namely that small scales 
tells us more about galaxy formation than it does about In other terms, the 

centers of galaxies are special places, the only places where we can study dark matter 
under peculiar conditions. 

4. Conclusions 

In this paper, we studied the cusp/core problem by improving the ESIM of Williams 
et al. (2004). We took into account, simultaneously, the effects of ordered and random 
angular momentum, dynamical friction and adiabatic contraction. The improved SIM of 
the present paper, taking account the previous effects gives rise to haloes being char- 
acterized by a log-log density slope that changes gradually within the virial radius and 
slopes of the inner part of haloes flattening with decreasing mass. As in previous papers, 
AC steepens the density profiles. The density profiles of structure having masses smaller 
than 10 11 hT x Mq are well fitted by Burkert's profiles. The comparison of some of the 
rotational curves given by de Blok & Bosma (2002) with the rotational curves obtained 
by means of our model, gives a good agreement. In the case of clusters of galaxies the 
density profile evolution is similar to that observed on galactic scales with the difference 
that the final slope is steeper than in the dwarf galaxies case. However the total mass 
profile is still cuspy and evolves slightly. The behavior of the dark matter halo is in 
agreement with the analysis of Sand et al. (2002, 2004) who fitted the baryonic and dark 
matter profiles only in the very inner part of the cluster MS 2137-23 within ~ dOh^ 1 kpc 
by means of a generalized NFW profile. The behavior of total mass is in agreement with 
X-ray observations by Chandra and XMM (Buote 2003, 2004; Lewis et al. 2003) weak 
lensing (Dahle et al. 2003) and strong lensing (Bartelmann 2003), which are consistent 
with a cusp having a = 1 or larger. As also reported by Williams et al. (2004), numeri- 
cally generated halos must have lost a considerable amount of their angular momentum 
in the outer parts, roughly between 0.1 and li?„, possibly through dynamical friction or 
other mechanisms. A further future model taking into account the interaction between 
DM and baryons would help to better understand the problem. 
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(b) log (R/RJ 

Fig. 1. Panel (a). Dark matter haloes generated with the model of Section 2. In this case 
we do not take into account baryon collapse. The upper dashed-line represents the results 
of the Aquarius Experiment for a halo of M ~ 1O 12 /i _1 M , the dot-dashed line represents 
the NFW profile for c = 10, calculated by means of Eqs. . (|4"7|) . (|4"5jl expressing the 
scaling radius r s in terms of the virial radius. The short-dashed line, dashed line and 
solid line represent, respectively, the density profile obtained by means of the model of 
the present paper for M ~ 1Q 12 /i _1 M q , M ~ 1O 1o /j _1 M , and M ~ lO 8 /!"^© without 
inclusion of baryonic infall. Panel (b) same of panel (a) but taking account of baryonic 
infall. 
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Fig. 2. Dark matter haloes generated with the model of Section 2. The dotted-dashed 
line is the NFW profile, the short-dashed line is the profile for the 10 12 /i _1 Mq halo 
obtained by means of our model. The dashed line represents the density profile obtained 
from the 1O 12 /i -1 M0 halo by reducing the magnitude of h, j and \i as described in the 
text. 
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Fig. 3. Comparison of the rotation curves obtained with the model in Section 2 (solid 
lines) with the rotation curves of four LSB galaxies studied by de Blok & Bosma (2002). 
The dashed line represents the fit with NFW model (see Section 3 for details). 
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Fig. 4. Panel (a) Density profile evolution of a 10 14 /i _1 Mo halo. The solid line represents 
the NFW initial profile at z = 3. Notice how the baryonic infall at redshift 2 < z < 3 
steepens the NFW cusp. The profile at z — 2, z — 1.5, z = 1 and z = is represented 
by the dot-short-dashed line, dotted line, short-dashed line, and long-dashed line respec- 
tively. Panel (b). Density profile evolution of a 10 14 ft- _1 M Q halo for the total mass. The 
solid line represents the NFW initial profile at z = 3. The profile at z — is represented 
by the dotted line. 



